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Nonequilibrium dynamical mean-field theory (DMFT) is developed for the case of the charge- 
density-wave ordered phase. We consider the spinless Falicov-Kimball model which can be solved 
exactly. This strongly correlated system is then placed in an uniform external dc electric field. 

We present a complete derivation for nonequilibrium dynamical mean-field theory Green’s functions 
defined on the Keldysh-Schwinger time contour. We also discuss numerical issues involved in solving 
the coupled equations. 
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I. INTRODUCTION 


There are a number of strongly correlated materials 
that have charge-density-wave (CDW) behavior. Static 
order occurs in the transition metal di- and trichal- 
chogenides, which display either quasi one dimensional 
(NbSea) or quasi two dimensional (TaSe 2 or TbTea) or¬ 
der—^. Three-dimensional charge-density-wave order is 
observed in Bai-^K^BiOa compounds— There is a long¬ 
standing question concerning the nature of the CDW or¬ 
der, namely, whether the order is driven electronically, 
with a lattice instability following the electronic instabil¬ 
ity, or vice versa. Indeed, in real materials the charge and 
lattice degrees of freedom are usually strongly coupled, 
but recent time-resolved core-level photoemission spec¬ 
troscopy experiments^ for some CDW materials indicate 
an electronically driven nature to ordering. This makes 
it reasonable to study the CDW phase for strongly cor¬ 
related electronic systems that do not include a coupling 
to the lattice. 

While most theoretical interests in strongly correlated 
systems have been concentrated on equilibrium behavior, 
recent experiments on pump-probe spectroscopy^ have 
caused an increase of attention to nonequilibrium dynam¬ 
ics. These experiments display a nonequilibrium melting 
of the CDW state, which is manifested by a filling of 
the gap in the photoemission spectrum, while the order 
parameter remains nonzero. This phenomenon has been 
examined with an exactly solvable model^ for a system 
starting at zero temperature. 

We use the Falicov-Kimball model in our analysis be¬ 
cause it is one of the simplest models^ which possesses 
static charge-density-wave ordering and has an exact so¬ 
lution within DMFT— (for a review see Ref—). Themany- 
body formalism for nonequilibrium dynamical mean-field 
theory is straightforward to develop within the Kadanoff- 
Baym-Keldysh formalism*2«— Since the many-body per¬ 
turbation theory diagrams are topologically identical for 
both equilibrium and nonequilibrium perturbation theo¬ 


ries^, the perturbative analysis of Metzner— guarantees 
that the nonequilibrium self-energy remains local. The 
basic structure of the iterative approach to solving the 
DMFT equations^ continues to hold. Detailed develop¬ 
ment of the nonequilibrium DMFT approach has been 
done for the case of the uniform phase of the Falicov- 
Kimball modeli— Here we generalize this method to the 
case of the CDW ordered phase. 


II. STATIC ORDER AND THE HAMILTONIAN 


In order to describe the CDW ordered state, one has 
to rewrite the Hamiltonian assuming the existence of the 
charge modulation. This can be done in two ways: by 
introducing two sublattices “A” and “B” in real space or 
by the nesting of the Brillouin zone (BZ) at the modu¬ 
lation vector Q = (n, 7 r,...) in a reciprocal space. The 
modulation vector Q defines the sublattices by 

e ‘ 9 R '=(-l, Tet (1) 


In the ordered phase, due to nesting of the Fermi sur¬ 
face, the BZ is reduced and instead of the annihilation 
(creation) operators with momentum k defined in the ini¬ 
tial BZ by Ck = 2 e* kRi Ci, one has to introduce two 
fermionic operators in momentum space in the reduced 
zone (k £ rBZ) as cm = Ck and 62 k = Ck+Q- 

Now, one can write down the relations between anni¬ 
hilation (creation) operators defined on the sublattices 
(A, B) and in the rBZ (1, 2) 
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or in a matrix form as follows: 
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III. REAL TIME GREEN’S FUNCTION IN THE 
CDW PHASE 

The key object of our interest is the time-dependent 
Green’s function that is defined on the Keldysh-Schwin¬ 
ger contou r 10 : 11 ’ 15 


We use both sublattice ( A,B ) and rBZ (1,2) bases. 
Any two-operator-product-type quantity, e.g. the one- 
particle Green’s function, can be defined with the ad¬ 
ditional sublattice indices <D( k) = ||O a ,/ 3 (k)|| (cq/3 = 

A,B), or with the rBZ indices 0( k) = ||O mi „(k)|| 
(■ m,n = 1,2), and these representations are connected 
by the aforementioned unitary transformation 


Gl(t,t') = -i(T c c k (t)cUt')). (8) 


We start from a formal solution of Dyson’s equation for 
the lattice Green’s function: 
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0(k) = Ud(k)U~ 1 . (4) 

We work in units where h = c = e = a= 1. 

The time-dependent Hamiltonian of the spinless Fali- 
cov-Kimball model on a bipartite lattice has the form 

U(t) = tifificlaCjP- ( 5 ) 
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We consider the case when charged fermions interact 
with an external electric field which is spatially uni¬ 
form. This allows us to describe the electric field via 
a time-dependent vector potential in the Coulomb gauge 
as E(i) = —dA(t)/dt. Interaction with this external field 
results in a Peierls’ substitution to the kinetic term of the 
Hamiltonian. 

In the sublattice representation (A, 13), the local part 
of the Hamiltonian is diagonal and the non-local kinetic 
one is off-diagonal in the case of the nearest-neighbor 
hopping. In the rBZ representation (1, 2) it is vice versa. 

The Fourier transformation to momentum space gives 
the time-dependent kinetic term in the form 
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where an extended band energy e(k — A(£))^ in matrix 
form in the (A, B) basis becomes 

e(k-A(£))= (7) 

0 e(k) cos(.A(£))+e(k) sin(A(t)) 

e(k) cos(.A(£))+e(k) sin0 

In the rBZ representation (1,2), an extended band 
energy matrix is diagonal. Here, we introduced the 
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In the two-sublattice representation (A,B), the self¬ 
energy E c {t,t') is diagonal and the noninteracting Gre¬ 
en’s function G('" orl '(f, £') is non-diagonal [because of ex¬ 
tended band energy e(k — A(£)) in Ecp 0], But in 
the rBZ representation (1,2) G c e ’2° n [t,t') becomes diag¬ 
onal and its analytical expression is known from the uni¬ 
form solution^. We apply the unitary transformation in 
Eq. 0 to the self-energy E c (f, t') and we find the solu¬ 
tion for the lattice Green’s function in Eq. 0 in explicit 
matrix form in the rBZ basis (1,2) as follows: 
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The Green’s functions and self-energies defined on the 
Keldysh-Schwinger contour are continuous matrix oper¬ 
ators of two time variables. We discretize the contour 
with several different grids and then find the continuous 
matrix limit using Lagrange’s interpolation formula. To 
find the inverse matrix in Eq. to, where its components 
are matrices in time variables, we need to apply the block 
matrix pseudo-inverse formula. 

We have performed a transition onto the BZ basis in 
order to find the lattice Green’s function in terms of the 
noninteracting Green’s function. Further, we construct 
the other DMFT equations to make the system of equa¬ 
tions self-consistent. In our case of a CDW phase, it is 
more convenient and clear how to do this in the sub¬ 
lattice basis (A,U). Hence, at this point we apply an 
inverse transformation from Eq. © onto the sublattice 
basis and write down expressions for the components of 
the lattice Green’s function in the (A, B) basis in terms 
of its components in the rBZ basis. 

The next step in the DMFT approach is to find the 
local Green’s function, which is further mapped onto a 
single impurity problem. Since we are considering a two- 
sublattice system, we have two local Green’s functions 
for the A and B sublattices, respectively. 

To calculate the local Green’s functions, we need to 
sum over the reduced BZ all (k — A(f))-dependent func¬ 
tions. Because of presence of the time-dependent vector 
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FIG. 1. Quadratic extrapolation of the current for U = 1.0 
with E = 1 at temperature T = 0.033 (An/ = 0.495). Differ¬ 
ent curves correspond to different At values. 


FIG. 2. Quadratic extrapolation of for U = 1.0 with E = 1 
at temperature T = 0.033 (An/ = 0.495). Different curves 
correspond to different At. 


potential, one has to replace the summation over the BZ 
by a double integration over the energies 

GfocM = = f de f *‘KM)G', e -(M') 

k 
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with a joint density of states which is a Gaussian in 
each variable for the hypercubic lattice^ and is given 
by p(e, e) = exp(—e 2 - e 2 )/-rr. 

Finally, we find the components of the local Green’s 
function for each sublattice as follows: 
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where we introduced new quantities K , A, and E which 
satisfy 
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Next, we need to map the local lattice Green’s function 
onto the impurity Green’s function. Employing Dyson’s 
equation, we introduce an effective medium 


Gf oc (t,t') = [(Gg)- 1 - E c ]"(14) 


The effective medium G^{t,t') is diagonal in the sublat¬ 
tice representation and its components are equal to 


G c 0 ’ a (t,t') = [(GZ)- 1 + (« = A,B) (15) 


On the other hand, it can be found from the Dyson’s 
equation that defines an effective dynamical mean field 
A c (t,t'). Its components for each sublattice equal to: 

G c 0 ’ a (t,t') = \{idt c + p d )5 c (tj) - A c '“(i,0] _1 - (16) 

We next extract the dynamical mean fields for each sub¬ 
lattice 

A.“■“(*, 0 = (W t c + p d )6 c (t,t') - (G c 0 ’ a r\t,t') (17) 

= (id t c + p d )S c (t,t') - (GZ c r\t lt ') - E C ’“(M'), 

which are the effective fields for the nonequilibrium 
single-impurity problems. 

Now, we close the system of DMFT equations with the 
solution of the impurity problem, which is known for the 
Falicov-Kimball model: 

GZ P M = [(1 -nJ)GZ+nJGZ}(tZ), where 

GlZt,t') = [(l-GZU)]- 1 GZ](t,t l ). (18) 

The difference between the A and B sublattices is de¬ 
fined by the order parameter An/, which is equal to the 
difference of the /-particle occupations at different sub¬ 
lattices (2An/ = n/ — nf). In the CDW phase, the 
total concentration of localized electrons is fixed, and 
An/ is defined from the initial equilibrium condition. In 
nonequilibrium, the order parameter remains the same 
as in equilibrium because the /-particles of the Falicov- 
Kimball model do not interact with the external field. 

We calculate the current in the Hamiltonian gauge by 
evaluating the operator average 

<j(*)> = + 0{t)F,t)G^(t,t), (19) 
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where the lesser Green’s function G^(t,t) is extracted 
from contour ordered Green’s function, and the velocity 
component is Uj(k) = lim t* sin(k/)/%/d. 
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The results of the numerical calculation depend stro¬ 
ngly on the discretization of the time interval and we 
need to take the limit of At —> 0. Smaller At results in 
larger matrices so we exploit an extrapolation procedure 
to get more accurate results. In Fig. [TJ we present the 
results for the current at U = 1.0 that corresponds to 
a metallic phase at temperature T = 0.033 and order 
parameter A rif = 0.495. Different curves correspond to 
different values of At and we use a quadratic Lagrange’s 
interpolation formula to extrapolate the result to At = 0. 
These show the correct zero current before the field is 
turned on at t = 0 (we choose E = 1) and their longer 
time behavior agrees with zero temperature calculations 
for the Falicov-Kimball model as well&. Similarly, in Fig. 
0 we show the results for concentration of d-electrons on 
the A-sublattice. Here, extrapolation shows an even more 
significant shift of the nonzero At results and its result 


is in a good agreement with the T = 0 calculations^. 
Other results for more values of the parameters will be 
presented elsewhere. 
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